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ABSTRACT 

The key features of the MATPHOT algorithm for precise and accurate stellar photometry and 
astrometry using discrete Point Spread Functions are described. A discrete Point Spread Func- 
tion (PSF) is a sampled version of a continuous PSF which describes the two-dimensional 
probability distribution of photons from a point source (star) just above the detector. The 
shape information about the photon scattering pattern of a discrete PSF is typically encoded 
using a numerical table (matrix) or a FITS image file. Discrete PSFs are shifted within an 
observational model using a 21 -pixel-wide damped sine function and position partial deriva- 
tives are computed using a five-point numerical differentiation formula. Precise and accurate 
stellar photometry and astrometry is achieved with undersampled CCD observations by using 
supersampled discrete PSFs that are sampled 2, 3, or more times more finely than the observa- 
tional data. The precision and accuracy of the MATPHOT algorithm is demonstrated by using 
the C-language mpd code to analyze simulated CCD stellar observations; measured perfor- 
mance is compared with a theoretical performance model. Detailed analysis of simulated Next 
Generation Space Telescope observations demonstrate that millipixel relative astrometry and 
millimag photometric precision is achievable with complicated space-based discrete PSFs. 

Key words: techniques: image processing, photometric — astrometry — instrumentation: 
detectors — methods: analytical, data analysis, numerical, statistical 



1 INTRODUCTION 

A Point Spread Function (PSF) is a continuous two-dimensional 
probability-distribution function which describes the scattering pat- 
tern of photons from a point source (star). 

Encoding a PSF as a continuous mathematical function works 
well for many ground-based astronomical observations due to the 
significant blurring caused by turbulence in the Earth's atmosphere 
and dome/telescope seeing. Ground-based PSFs are typically char- 
acterized by having a lot of the power in their spatial-frequency 
distributions at low spatial frequencies. 

Space-based PSFs frequently have significant amounts of 
power at higher spatial frequencies due to the lack of blurring 
caused by atmospheric turbulence. Adaptive optics can produce 
PSFs with characteristics found in both uncorrected ground-based 
PSFs and space-based PSFs: low-spatial-frequency features (e.g., 
broad halos) are frequently combined with high-spatial-frequency 
features (e.g., due to segmented mirrors). 

Some PSF-fitting stellar photometric reduction programs de- 
scribe the PSF as a combination of continuous mathematical func- 
tions and a residual matrix which contains the difference between 
the mathematical model of the PSF and an observed ("true") PSF. 
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This artificial breaking of the PSF into analytical and discrete com- 
ponents is not without mathematical risk. Such residuals can have 
small features which are described with higher spatial frequencies 
than are present in the actual observational data — a problem that 
can usually be mitigated by sampling residuals at higher spatial res- 
olutions than the observational data. 

What if we dispose of the use of continuous mathematical 
functions to model any part of the PSF and just use a matrix to 
describe all of the PSF? Is precise and accurate stellar photome- 
try and astrometry possible using matrix PSFs with oversampled 
stellar image data? If that is possible, then what extra information, 
if any, is required in order to do precision photometric reductions 
with matrix PSFs on undersampled data? 

This article describes how precise and accurate stellar pho- 
tometry may be obtained using PSFs encoded as a matrix. The fol- 
lowing section derives the theoretical performance limits of PSF- 
fitting stellar photometry and astrometry. Some of the key features 
of the MATPHOT algorithm are presented in §3. A demonstra- 
tion computer program, called mpd, based on the current imple- 
mentation of the MATPHOT algorithm, is described in §4. Sim- 
ulated CCD (charge-coupled device) stellar observations are ana- 
lyzed with mpd in §5 and the performance of the MATPHOT al- 
gorithm is compared with theoretical expectations. Concluding re- 
marks are given in §6. An appendix explains box-and-whisker plots 
which are used extensively in this article. 
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2 THEORETICAL PERFORMANCE LIMITS 
2.1 Point Response Functions 

A Point Response Function, <]/, is the convolution of a Point Spread 
Function, <f>, and a Detector Response Function, A : 



* = 0* A . 



(1) 



The PSF describes the two-dimensional distribution of photons 
from a star just above the detector. Although stellar photons are 
distributed as a point source above the Earth's atmosphere, a stellar 
image becomes a two-dimensional distribution as the stellar pho- 
tons are scattered by atmospheric turbulence. The blurred stellar 
image is then further degraded by passage of the stellar photons 
through the combined telescope and camera optical elements (such 
as mirrors, lenses, apertures, etc.). The PSF is the convolution of 
all these blurring effects on the original point-source stellar image. 
The two-dimensional discrete (sampled) Detector Response Func- 
tion (DRF) describes how the detector electronics convert stellar 
photons (7) to electrons (e~) — including such effects as the dif- 
fusion of electrons within the detector substrate or the reflection 
(absorption) of photons on (in) the gate structures of the detector 
electronics. 

The PSF is a two-dimensional probability-distribution func- 
tion describing the scattering pattern of a photon. The volume in- 
tegral of the PSF is one: Vpsf = 1; photons, after all, have to be 
scattered somewhere. It is important to note that since the angular 
extent of a PSF can be quite large, the volume integral the PSF over 
any given observation is frequently less than one due to the limited 
spatial coverage of the observation. 

The volume integral of a PRF is, by definition, one or less: 
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where a value that is less than one indicates a loss of stellar photons 
during the detection/conversion process within the detector. While 
the quantum efficiency (QE) variations within a single detector are 
generally not a major problem with state-of-the-art charge-coupled 
devices, intrapixel QE variations can be significant with some near- 
infrared detector technologies currently being used in astronomical 
cameras (e.g., Lauer 1999, Hook & Fruchter 2000). 

A perfect DRF gives a PRF that is a sampled version of the 

PSF: 

$i= / cj>(x,y)dxdy , (3) 

Jxi-0.5 Jyi-0.5 

where i th pixel of the PRF located at (Xi,yt) is the volume inte- 
gral of the PSF over the area of the i pixel. The actual limits of 
the above volume integral reflect the appropriate mapping transfor- 
mation of the x and y coordinates onto the CCD pixel coordinate 
system. 

The sharpness of a PRF is defined as the volume integral of 
the square of the normalized PRF: 
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Physically, sharpness is a shape parameter which describes the 
"pointiness" of a PRF; sharpness values range from a maximum 
of one (all of the stellar flux is found within a single pixel) to a 
minimum of zero (a flat stellar image). For example, cameras that 



are out of focus have broad PSFs with sharpness values near zero. 
A normalized Gaussian PSF with a standard deviation of S pixels, 
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that has been oversampled with a perfect DRF will have a 
sharpness value of 



g (x,y;X,y,S) dxdy 



1 



(6) 



A critically-sampled normalized Gaussian PRF has a sharpness of 
l/(4n) and any PRF with a sharpness value greater than that value 
(~0.0796) can be described as being undersampled. Diffraction 
limited optics, theoretically, give sharpness values that decrease 
(i.e., PSFs become flatter) with increasing photon wavelength - for 
a fixed pixel (detector) size. With real astronomical cameras, the 
value of sharpness frequently depends on where the center of a star 
is located within the central pixel of the stellar image. For exam- 
ple, the Hubble Space Telescope (HST) WFPC2 Planetary Camera 
PRF at a wavelength of 200 nm has an observed sharpness value 
of 0.084 if the PRF is centered in the middle of a PC p ixel or 0.063 
if the PRF is centered on a pixel corner (Table 6.5 of lBiretta et alJ 
2001); at 600 nm the observed sharpness values range from 0.066 
(pixel-centered) to 0.054 (corner-centered). The Wide-Field Cam- 
eras of the HST WFPC2 instrument have pixels which are approx- 
imately half the angular resolution of the PC camera pixels; stel- 
lar images on the WF cameras are undersampled and the observed 
range of WF camera sharpness values are 0.102-0.120 at 200 nm 
and 0.098-0.128 at 600 nm. 

The effective-background area, (3, of a PRF is defined as the 
reciprocal of the volume integral of the square of the PRF: 



[3 = 
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Alternatively, the effective-background area (a.k.a. equivalent- 
noise area or effective solid angle ) of a PRF is equal to the recip- 
rocal of the product of its sharpness and the square of its volume: 
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The effective-background area of a normalized Gaussian PRF is 
47TiS 2 px 2 , where S is the standard deviation in pixels; a critically- 
sampled normalized Gau ssian PRF has an effective-background 
area of 47r w 12.57 px. iKina 11983) notes that numerical inte- 
gration of a realistic ground-based stellar profile gives an effective- 
background area of 30.8 <S 2 instead of the value of 47rS 2 for a 
normalized Gaussian profile. 

2.2 Basic Least-Squares Fitting Theory 

Consider a CCD observation of two overlapping stellar images. As- 
suming that we already know the PSF and the DRF of the observa- 
tion, a simple model of the observation will have seven parameters: 
two stellar intensities 1 {£1,82) in electrons, four coordinate val- 

1 Stellar intensity is denned to be the total number of electrons from a sin- 
gle star scaled to a PRF volume integral of one. The observed stellar inten- 
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ues, giving the stellar positions (Xi, 3^i, <-f 2, 3^2) in pixels, and B 
which is the observed background sky level 2 in electrons (which is 
assumed to be the same for both stars). These observational param- 
eters are not independent for overlapping stars in the presence of 
photon and CCD readout noise. The conservation of electron flux 
will require that if £\ increases then £2 must decrease and vice 
versa for a given value of B. The most accurate photometry possi- 
ble is obtained when these dependent parameters are fitted simul- 
taneously. Any reasonable model of two overlapping stellar images 
will be a nonlinear function when the positions and intensities are 
to be determined simultaneously. The technique of nonlinear least- 
squares fitting was developed to provide for the simultaneous de- 
termination of dependent or independent parameters of nonlinear 
model functions. 

Assume that we have a calibrated CCD observation with N 
pixels and that Zi is the number of electrons in the i th pixel which 
is located at the position of (xi,yi) and has a measurement error of 
o i electrons. Let m(x,y;pi, . . . ,Pm) be an observational model 
of the CCD electron pixel values that has two coordinates (x, y) 
and M parameters. For notational convenience, let the vector r; 
represent the coordinates (xi,yi) of the i th pixel and the vector 
p represent all the model parameters [p = (pi, . . . ,Pm) ]. The 
observational model of the i th pixel can thus be compactly written 
as follows: m 4 = m(r 4 ; p). 

The measure of the goodness of fit between the data and the 
model, called chi square, is defined as 



x 2 (p) 
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(9) 



The theory of least-squares minimization states that the optimum 
value of the parameter vector p is obtained when x 2 (p) ls mini- 
mized with respect to each parameter simultaneously. If p is the 
optimal parameter vector, then X 2 (Po) ' s me absolute minimum of 
the A/-dimensional manifold X 2 (p)- 

For some small correction parameter vector S one can approx- 
imate x 2 (p + 8) by its Taylor series expansion: 
00 
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where 
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is th e jk th element of the MxM H essian matrix H of x 2 (p) 
[e.g.,|Arfken 1 1970); P ress et alJ ll986)l. The approximation for the 
calculation of the Hessian matrix elements is frequently used when- 
ever the computation of the second partial derivative is numerically 
unstable. If x 2 (p + 8) is a local minimum of x 2 manifold, then it 
can be shown that 



H-S = -V X 2 (P) 



(12) 



sity (= EN) is, by definition, is always less than or equal to the measured 
stellar intensity {= £) . 

2 The observed background sky level (in electrons) is the product of true 
background sky level (in photons) and the average PRF volume across a 
pixel: B = Strue(V) . 



By solving this equation for the correction vector S , one can deter- 
mine a better parameter vector as follows: p' = p + S . When the 
parameter vector (p) is redefined to be the better parameter (p'), 
the Hessian matrix and the gradient of x 2 (p) can then be recalcu- 
lated to determine a new correction vector (S). This process repeats 
until the correction vector is sufficiently small - generally when the 
difference between solutions is no longer statistically significant. If 
the fitting process has not failed, then the optimal parameter vector 
(p ) should be very close to the true parameter vector. 

Once the optimal parameter vector has been determined, the 
covariance matrix C may then be calculated by inverting the Hes- 
sian matrix H computed with the optimal parameter vector. The 
standard errors (one standard deviation) of the fitted parameters can 
be estimated as follows: 
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where aj is the standard error associated with the j parameter 
(Pj). Usage of equation 1 1 3i for error estimates is based on the crit- 
ical assumption that fitted model parameters are independent (indi- 
cated by negligibly small off-diagonal elements of the covariance 
matrix). It is important to note that whenever this critical assump- 
tion is violated, the results produced by least-squares fitting may 
not be statistically reliable, which is to say, they may no longer be 
physically meaningful. 



2.3 Photometry 

The theoretical photometric performance limits for PSF-fitting 
CCD stellar photometry can be derived using a simple observa- 
tional model consisting of a PRF and a constant sky level. 



2.3.1 Observational Model 

Consider a CCD observation of single isolated star on a flat sky 
background. Assuming one already knows the PRF of the observa- 
tion at the location of the star, a simple model of the observation 
would have just two parameters: the stellar intensity (£) in elec- 
trons, and the observed background sky level (B) in electrons. The 
observational model for the i th pixel would be 



(14) 



where V is the volume integral of the PRF and $j is the value of 
the i th pixel of the normalized PRF ( = *i/V ). 



2.3.2 Bright Star Limit 

In the case of bright stars, most of the electrons found in the i 
pixel of the observation will come from the star and not the sky: 

rm « £V*» . (15) 

The actual number of electrons found in the i th pixel will be de- 
scribed by a Poisson distribution with a mean and variance of mj. 
The measurement error (one standard deviation) for the i th pixel 
would thus be 



Oi = 



(16) 
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All other noise sources (due to, for example, the observed back- 
ground sky, instrumental readout noise, flat-field calibrations er- 
rors, etc.) are assumed, in this special case, to be negligibly small. 

The variance of the stellar intensity measurement error of 
bright stars can be estimated using equations i 1 3i . 1 141 . and 1161 : 



^bright 
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dx dy 
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as expected from photon statistics. 

A bright isolated star with an intensity of 10 6 photons im- 
aged with a perfect CCD detector would have a stellar image 
with 10 6 e~ (= £) and a stellar intensity measurement error 
of ag ss y/E /(V = 1) = 10 3 e~. The same star imaged with 
an inefficient CCD detector with a quantum efficiency of 25% 
(V = 1/4) would have a stellar image with ~250, 000 e - which 
would have a Poisson noise error of ~500 e" . The measured stel- 
lar intensity is £ ~ 10 6 e~ with an rms measurement error of 
erg ~ \ ft~fV = 2000 e - which is two times larger than it would 
be with a perfect detector and four times larger than the Poisson 
noise error of the observed number of electrons. 

Solving for measured stellar intensity (= £) instead of the 
observed stellar intensity (= £ V) enables the creation of stellar 
photometric reduction programs capable of dealing with intrapixel 
QE variations through the accurate modeling of the image forma- 
tion process within the detector. While it is certainly convenient to 
assume that one's detector has negligible intrapixel QE variation, in 
the real world even NASA-grade CCD detectors, like those found 
in the HST WFPC2 instrument, can have peak-to-peak intrapixel 
sensitivity vari ations that a re greater than 0.02 mag (>2%) (see 
Figs. 5 and 6 of lLaueJl999l) . 



2.3.3 Faint Star Limit 

Most of the electrons found in the i th pixel of an observation of a 
faint isolated star on a flat sky background will come from the sky 
and not from the star. In that case, the measurement error associ- 
ated with i th pixel is approximately the effective-background noise 
level: 



where 
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B is the constant observed background sky level which is assumed 
to be a Poisson distribution with a mean of B electrons, and ctron 
is the rms readout noise. 

The variance of the stellar intensity measurement error of faint 
stars can be estimated using equations 1131 . i 141 . 1 1 81 . 1 1 9i . I20i . 



and©: 
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where (3 is the effective-background area of the PRE Equation 1221 
agrees with equation (9) of Kina 1 1983) for a perfect (V= 1) noise- 
less (ctron = e~ ) detector. 

An important additional noise source for the photometry of 
faint stars is the systematic error due to the uncertainty of the mea- 
surement of the background. If the sky background is assumed to 
be flat, then the rms measurement error of the constant sky back- 
ground can be estimated using equations 1131 . 1 141 . d 1 8i . 1 1 91 . and 
03: 
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Given a CCD observation with no readout noise, equation 1241 re- 
duces to the value of erg = yjB/N expected from simple sam- 
pling statistics. 

The portion of the rms stellar intensity measurement error that 
is cause d by the erro r in the determination of the local sky level 
is <rg (3 (Irwin 1985). While this error is frequently negligible for 
bright stars, it is generally significant for faint stars. Including the 
uncertainty in the determination of the constant observed back- 
ground sky level thus gives a more realistic estimate for the rms 
stellar intensity measurement error for faint stars: 
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Precise and accurate stellar photometry of faint stars requires an 
excellent determination of the observed background sky which in 
turn requires accurate background sky models. Given a valid back- 
ground sky model, small apertures will be more sensitive to back- 
ground sky measurement errors than large apertures. 



2.3.4 Photometric Performance Model 

A realistic photometric performance model for CCD PSF-fitting 
photometry can be created by combining the bright and faint star 
limits developed above. The theoretical upper limit for the photo- 
metric signal-to-noise ratio (S/N) of CCD PSF-fitting photometric 
algorithms is as follows: 
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T £: bright ^ G £: faint 
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These approximations assume, for the sake of simplicity, that any 
noise contribution due to dark current and quantization noise is neg- 
ligible. While these additional noise sources can be added to cre- 
ate an even more realistic performance model for stellar photom- 
etry, the assumption of low dark current and minimal quantization 
noise is realistic for state-of-the-art astronomical-grade CCD im- 
agers. The resulting photometric error is approximately 



Amag 



1.0857 
S/N 



(29) 



where the constant 1.0857 is an approximation for Pogson's ratio 
as5/ln(100)=2.51og(e) JPogsoi]ll85d) . 



2.3.5 Cramer-Rao Lower Bound 

The Cramer-Rao Lower Bound (CRLB) is the lower bound on the 
variance of any unbiased estimator. Since it is physically impossi- 
ble to find an unbiased estimator that beats the CRLB, the CRLB 
provides a performance benchmark against which any unbiased es- 
timator can be compared. 

The Cramer-Rao Lower Bound for stellar photometry of 
a single isolated star imaged by a two-dimensional photon- 
counting detector has been derived several ti mes in the as- 
troph y sical l i teratur e (see , e.g., Appe ndix A of IPerrvman et alJ 
I1989L llrwirJ Il985l and iKind Il983l) . The generalization for 
a crow ded fiel d with overlapping stellar images is given in 
Ijakobsen. Greenfield. & Je drzeiewskJ ll992l) . 

The Cramer-Rao Lower Bound for the bright star limit of stel- 
lar photometry of a single isolated star is 



J £:bright-CRLB 



= £ 



(30) 



which is equation I17i with a perfect detector. The CRLB for the 
faint star limit of stellar photometry of a single isolated star is 
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a £ 



f3B 



(31) 



faint-CRLB 

which is equation 1261 with a noiseless detector and a negligible 
background measurement error (N^oo). 

The photometric performance model has bright and faint star 
limits which are the same, respectively, as the bright and faint star 
Cramer-Rao Lower Bounds for stellar photometry of a single iso- 
lated star on a flat sky background imaged with a perfect noiseless 
detector. 



2.4 Astrometry 

The theoretical astrometric limits for PSF-fitting CCD stellar pho- 
tometry can be derived using a simple observational model consist- 
ing of a Gaussian PRF and a constant sky level. 



2.4.1 Observational Model 

Consider a CCD observation of single isolated star on a flat sky 
background. A Gaussian is a good model for the PSF of a ground- 
based CCD observation since the central core of a g round-based 
stellar profile is approximately Gaussian (King 1971). In this case 
the PSF would have three parameters: two coordinate values giv- 
ing the location (X,y) of the star on the CCD and the standard 
deviation of the Gaussian (5) in pixels [see equation 

An imperfect but uniformly flat DRF (V < 1) gives a value 
for the i th pixel of the PRF located at (a;,-, yi) which is equal to 
the product of the volume of the PRF and the value of the volume 
integral of the PSF over the area of the i pixel: 



G ; sV 
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g(x,y;X,y,S) dxdy . 



(32) 



The actual limits of the above volume integral reflect the appro- 
priate mapping transformation of the x and y coordinates onto the 
CCD pixel coordinate system. 

If the PRF has been oversampled, the value of the i pixel of 
the PRF is approximately equal to the product of the volume of the 
PRF and the value of the PSF at the center of the i th pixel: 



G^V gl 
where 

g 1 = &(xi,yi;X,y,S) 



(33) 



(34) 



A simple model of the observation will require two additional 
parameters: the stellar intensity (£) in electrons and the observed 



background sky level (B) in electrons. The i 
vational model would be 

rm = B + £VG Z , 



pixel of the obser- 



(35) 



where V is the volume integral of the PRF and Gi is the value of 
the i th pixel of the normalized PRF (Gi = Gi /V ~ gi). 



2.4.2 Bright Star Limit 

In the case of bright stars, most of the electrons found in the i 
pixel of the observation will come from the star and not the sky: 



£VGi 



(36) 



The actual number of electrons found in the i pixel will be de- 
scribed by a Poisson distribution with a mean and variance of rrii. 



The measurement error (one standard deviation) for the i 
would thus be 

Ci — J~mi 



pixel 



(37) 



All other noise sources (e.g., due to the observed background sky, 
instrumental readout noise, flat-field calibrations errors, etc.) are 
assumed to be negligibly small. 

The variance of the stellar X position measurement error of a 
bright isolated oversampled star can be estimated using equations 
O- 03, OIK and J5}: 



J X: bright 
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(39) 



is the critical-sampling scale length of the PRF 3 in pixel units (px), 
which, unlike S, is defined for all PRFs. By definition, the critical- 
sampling scale length of a critically-sampled PRF imaged with a 
perfect detector is one pixel; C > 1 indicates that the PRF is over- 
sampled, while C < 1 indicates that the PRF is undersampled. 

In the special case of a critically-sampled bright star imaged 
with a perfect detector, one finds that the astrometric performance 
limit (in pixel units) is equal to the reciprocal of photometric error 
performance limit: 



bright 



V£ a £: bright 



2.4.3 Faint Star Limit 

Let us again assume that the noise contribution from the star is 
negligibly small and that the variance of the measurement error of 
the i th pixel can be replaced with an average constant rms value. 
The variance of the stellar X position measurement error of a faint 
isolated oversampled star can be estimated using equations 1 1 3i . 
{35j, HD, 033. HD, and : 
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(40) 
(41) 



3 From the definition of the effective-background area of an oversampled 
Gaussian PRF with V < 1, @ G = 4tt5 2 /V 2 , one sees that critical- 
sampling scale length has been designed to be a proxy for S for any PRF. 



2.4.4 Astrometric Performance Model 

A realistic performance model for CCD PSF-fitting astrometry can 
be created by combining the bright and faint star limits developed 
above. The expected lower limit of the rms measurement error for 
the stellar X position for a single isolated star on a flat sky can be 
estimated as follows: 



a X 



1 1 
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1 + 8^ (B + <4on) 
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(43) 



The rms stellar y position measurement error is, by symmetry, the 
same as for X : 

a y ~ a X ■ ( 44 ) 

2.4.5 Photonic Limit and the Cramer-Rao Lower Bound 

The Cramer-Rao Lower Bound for stellar astrometry depends not 
only on the signal-to-noise ratio but also on the size and shape of 
the detector. For well-sampled data, the size and shape of the de- 
tector can be ignored and a CRLB can be found for a perfect noise- 
less detector with infinitely small pixels. This is called the photonic 
limit. 

The determination of the CRLB for astrometry becomes much 
more complicated with undersampled observations. Astrometric 
precision degrades when the size of the detector is comparable to 
the size of the stellar image - the quality of the position estimation 
is then dependent on the fraction of photons falling outside of the 
central pixel. The worst-case scenario for stellar astrometry occurs 
when all of the light from a star falls within a single pixel: all one 
knows for sure, in that unfortunate case, is that the star is located 
somewhere within the central (and only) pixel. 

The photonic limit (PL) for stellar astrometry of a bright well- 
sampled single isolated normalized Gaussian star is 

2 S 2 

a X: bright-PL = ~g 

(Irwin 1985). Using £ as a proxy for <S, one has the generalized 
form for any PSF: 



(45) 



2 £ 2 

a X: bright-PL ~ ~g ' 

which is equation J38> with a perfect detector. 

The photonic limit for stellar astrometry of a faint well 
sampled single isolated normalized Gaussian star is 



T X: faint-PL 



8nBS 4 
£ 2 



llrwinlll985t) . Using C as a proxy for 5, one has the generalized 
form for any PSF: 

2 8ttS£ 4 



a X: faint-PL ; 



£ 2 



(46) 



which is equation 1-4 1 i with a perfect noiseless detector. 

The astrometric performance model has bright and faint star 
limits which are the same, respectively, as the bright and faint 
star photonic astrometric limits, which are the Cramer-Rao Lower 
Bounds for stellar astrometry of a single isolated Gaussian star on 
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a flat sky background imaged with a perfect noiseless detector with 
infinitely small pixels. The Cramer-Rao Lower Bound for stellar 
astrometry of a single isolated Gaussian star on a flat sky back- 
ground imaged with a perfect noiseless CCD with square pixels 
<Winicldfl98a) quickly approaches the photonic limits with well- 
sampled observations; undersampled observations will have larger 
astrometric errors than predicted by the photonic limits. 



2.5 Relation between Astrometric and Photometric Errors 

2.5.1 Bright Star Limit 

Following Elngl <1983l) and llrwinl 1 19851) . I now compare the as- 
trometric error of bright isolated stars with their photometric error. 
The ratio of the astrometric error of a bright isolated star and the 
critical-sampling scale length of the PRF is equal to the ratio of the 
stellar intensity measurement error and the stellar intensity: 



ox 
C 



1L 
£ 



(47) 



For example, a bright isolated critically-sampled star with one mil- 
lion electrons imaged on a perfect detector (£ — 10 6 e _ , V = 1 , 
and C — 1 px ) would, theoretically, have a signal-to-noise ratio of 
a thousand (S/N = 1000), a stellar intensity measurement error of 
og = 1000 e~ , and an rms position error in x of one-thousandth 
of a pixel (p x = 0-001 px). Such astrometric accuracy may be 
difficult to achieve in practice under normal ground-based observ- 
ing conditions even with state-of-the-art astronomical-grade CCD 
cameras. 



2.5.2 Faint Star Limit 

The astrometric error of faint isolated stars is related to their pho- 
tometric error as follows: 



ox ^ ( o £ \ 
T ~ V £ ) T 



V2 



(48) 



For example, a faint isolated critically-sampled star imaged with a 
perfect detector with a 20.0% intensity measurement error and a 
negligible background measurement error (iV — ► oo) would, theo- 
retically, have an astrometric error of ~0.283 [~ (0.200) \pi ] px . 



2.5.3 Practical Lower Bound 

These results suggest the following practical lower bound for astro- 
metric errors with respect to photometric errors: 

X % photometry gives no better than X % astrometry 
with respect to the critical-sampling scale length (£.). 

For example, a star with one-percent stellar photometry will have 
no better than one-percent astrometry with respect to the critical- 
sampling scale length. If the star is critically sampled, then the as- 
trometric precision will be no better than 0.01 px. 

All of the above derivations are based on the assumption that 
that flat-field calibration errors are negligible. The relation between 
photometry and astrometry for bright isolated stars can fail with 
large flat-field calibration errors. 



3 DISCRETE POINT SPREAD FUNCTIONS 

A discrete Point Spread Function is a sampled version of a con- 
tinuous two-dimensional Point Spread Function. The shape infor- 



mation about the photon scattering pattern of a discrete PSF is 
typically encoded using a numerical table (matrix). An analyti- 
cal PSF has the shape information encoded with continuous two- 
dimensional mathematical functions. 

In order to do accurate stellar photometry and astrometry with 
discrete PSFs one needs to able to (1) accurately shift discrete PSFs 
to new positions within the observational model, and (2) compute 
the position partial derivatives of discrete PSFs. The next two sub- 
sections describe how these tasks may be accomplished using nu- 
merical analysis techniques. 



3.1 Moving Discrete PSFs 

Building a realistic observation model requires the placement of a 
star at the desired location within the model; this is done by de- 
termining the PRF at required location and then multiplying it by 
the stellar intensity. With PSFs encoded by mathematical functions, 
one just computes the PSF at the desired location in the observa- 
tional model. With discrete PSFs, one ideally takes a reference PSF 
(typically derived/computed for the center of a pixel) and shifts it 
to the desired location using a perfect two-dimensional interpola- 
tion function. But how is this done in practice? The sine function, 
sin(7ra;) / (irx), is, theoretically, a perfect two-dimensional interpo- 
lation function. Unfortunately, the sine function decays with 1/x 
and never actually reaches zero. One can use a windowed inter- 
polant in order to improve computational speed — but one must be 
cautious about aliasing effects caused by using a windowed func- 
tion. In the case of stellar photometry and astrometry, aliasing ef- 
fects will generally only be seen with bright stars since a large num- 
ber of photons are required in order to adequately sample the higher 
spatial frequencies of the PSF. 

The following 21 -pixel- wide damped sine function interpolant 
does an excellent job interpolating discrete PSFs: 



r 



J (*o) 



10 



i=-10 



^sin (ii(xi — xq)) 
n(xi — x ) 



exp 



Xo 



3.25 



(49) 



(Mishell 2002). Note that since the two-dimensional sine function 
is separable in x and y, this interpolant can be coded to be compu- 
tationally fast and efficient. This interpolant, from the ZODIAC C 
library written by Marc Buie of Lowell Observatory, was specifi- 
cally designed for use with 32-bit floating numbers. 

Aliasing problems due to critically-sampled or undersampled 
data may be overcome by using discrete PSFs that are super- 
sampled at 2, 3, or more times more finely than the observational 
data. In order to have a realistic observational model, once the su- 
persampled discrete PSF has been interpolated to the correct posi- 
tion, a new degraded (rebinned) version of the discrete PSF must be 
created which has the same spatial resolution as the observational 
data. 



3.2 Position Partial Derivatives of Discrete PSFs 

While the mathematics of determining the position partial deriva- 
tives of individual stars within the observational model with re- 
spect to the x and y direction vectors is the same regardless of 
how the shape information in a PSF is encoded, the implementa- 
tion methodology for the computation of position partial derivatives 
of discrete PSFs is very different than the one used for analytical 
PSFs. 
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The position partial derivatives of discrete PSFs can be de- 
termined using numerical differentiation techniques on the discrete 
PSF. 

It is a standard practice in numerical analysis to approximate 
the first, second, or higher, deriva tives of a tabulated functi on f(xi) 
with multi-point formulae. Abra mowitz & Stegunl 1196 4) give 18 
different multi-point formulae which can be used (with varying de- 
grees of accuracy) to approximate the first derivative of the tabu- 
lated function f(xi). The following five-point differentiation for- 
mula, 

/'(**) 

m ^ [/(au-a) ~ 8 f(xi-x) + 8 f(x l+1 ) - f(x i+2 )] , (50) 

(p. 91 4 of lA bramowi tz & Stegur]| 19641) works well with discrete 
PSFs iMiahell 2002). This approximation takes just 4 additions and 
3 multiplications which generally makes it considerably faster to 
compute than the traditional determination of the partial derivative 
of the volume integral of the PSF above a CCD pixel. 



4 THE MATPHOT ALGORITHM 

The concepts presented above outline the unique and fundamen- 
tal features of the MATPHOT algorithm for accurate and precise 
stellar photometry using discrete Point Spread Functions. 

While the key features of a CCD stellar photometric reduc- 
tion algorithm can be described in an article, the full implementa- 
tion of such an algorithm generally exists as a complex computer 
program consisting of many thousands of lines of computer code. 
Since good algorithms can be poorly implemented, it can be diffi- 
cult to differentiate between a poor algorithm and a poorly-coded 
implementation of a good algorithm. 

Confidence in a complex algorithm can be established by de- 
veloping an implementation of the algorithm which meets theoreti- 
cal performance expectations. The following subsection describes a 
real-world implementation of the MATPHOT algorithm that meets 
the theoretical performance expectations for accurate and precise 
stellar photometry and astrometry which were derived in §2. 

4.1 MPD: MatPhot Demonstrator 

I have written a C-language computer program, called mpd 4 , 
which is based on the current implementation of the MATPHOT al- 
gorithm for precise and accurate stellar photometry using discrete 
Point Spread Functions. The mpd code demonstrates the precision 
and accuracy of the MATPHOT algorithm by analyzing simulated 
CCD observation s based on user-prov ided discrete PSFs encoded 
as FITS images IWells. Greisen. & Hartei]|l98ll) . Discrete PSFs 
are shifted within the observational model using the 21 -pixel- wide 
damped sine interpolation function given in equation I49i . Position 
partial derivatives of discrete PSFs are computed using the five- 
point differentiation formula given in equation 1501 . Accurate and 
precise stellar photometry and astrometry of undersampled CCD 
observations can be obtained with the mpd code when it is pre- 
sented with supersampled discrete PSFs that are sampled 2, 3, or 
more times more finely than the observational data. The mpd code 
is based on a robust implementation of the Levenberg-Marquardt 

4 All source code and documentation for mpd and support soft- 
ware is freely available at the official MATPHOT website at NOAO: 
http://www.noao.edu/staft7mighell/matphot 



method of nonlinear least-squares minimization jLevenberall 19441 
lMarauarddfl963l also lMighelll ll989). When presented with simu- 
lated observations based on a Gaussian PSF with a known FWHM 
value °, the mpd code can analyze the observation two different 
ways: (1) the MATPHOT algorithm can be used wi th a discrete 
Gaussian PSF, or (2) analytical techniques lMighertfl989l.il 9991) 
can be used with an analytical Gaussian PSF. 



5 SIMULATED OBSERVATIONS 
5.1 Oversampled PSFs 

I now demonstrate that the theoretical performance limits of §2 pro- 
vide practical performance metrics for photometry and astrometry 
of CCD stellar observations that are analyzed with oversampled 
Gaussian Point Spread Functions. 

5.1.1 Analytical PSFs 

Twenty thousand oversampled CCD stellar observations were sim- 
ulated and analyzed using the mpd code. The CCD detector was 
assumed to be perfect (V = 1) with a CCD readout noise value 
of ctron = 3 e~ px -1 . Stars were simulated using an analytical 
Gaussian PSF with a FWHM = 3 px located near the center of 
60x60 pixels, the input stellar intensities ranged from —6 to —15 
mag 6 (251 ^ £truc ^ 10 6 e~), and a flat background was assumed 
with a value of B — 100 e - . Photon and readout noise was simu- 
lated, respectively, using Poisson and Gaussian random noise gen- 
erators and the resulting observed background sky measurement 
error was erg = 0.18 e~. The median effective-background area 
of the PRF of these observations was /3 = 21.44 px 2 . All the sim- 
ulated observations were analyzed with mpd using an analytical 
Gaussian PSF with FWHM = 3.0 px. 



The binned absolute photometric errors are shown as black 
box-and- whiskers plots (see AppendixlA"! on the top panel of Fig.Q 
The absolute photometric error of an observation is the absolute 
value of the difference between the measured (estimated) and true 
(actual) stellar magnitude: Amag = |mag — mag truc |. The four 
grey limits seen on the top panel of Fig. Q are theoretical predic- 
tions (derived from §2.3.4) for the median (50% cumulative frac- 
tion: grey solid curve), top hinge (75%: bottom of the grey band), 
top fence (~98.35%: top of band), and 5-cr outlier (~99.99997%: 
grey dashed curve) values. If the rms photometric error is called 
fmag, then the values of these theoretical limits are approximately 
equal to 0.674 amag, 1.151 ornag, 2.398 crmag, and 5.423 ornag, 
respectively. If the photometric performance model is correct and 
mpd has been coded correctly, then (1) the observed median val- 
ues (central bar in each box) should intersect the theoretical median 
value, (2) most of the top whiskers should be found inside the band, 
and (3) most of the outliers should be found above the top of the 
band and all of the outliers should found below the 5-cr outlier limit. 

5 The FWHM (Full-Width-at-Half-Maximum) value of a Gaussian is 
equal to 2^/ln(4) times the standard deviation, S, of the Gaussian: 
FWHM M 2.35482 5 [see equation 0], 

6 The MATPHOT magnitude system assumes that mag = 1 e — (electron) 
= 17 (photon) for a Point Response Function volume of one (V = 1). 
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Figure 1. The absolute photometric errors (top) and total astrometric er- 
rors (bottom) of 20, 000 simulated CCD stellar observations analyzed with 
mpd using an oversampled analytical Gaussian PSF with a FWH M of 3.0 
px(/3 ss 21.44 px 2 ; V = 1). 
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Figure 2. Relative stellar intensity errors (top) and relative X position errors 
(bottom) of the data set used in Fig. HI 



Comparing the absolute photometric errors of the 20,000 sim- 
ulated CCD observations with the grey theoretical limits, one sees 
that the photometric performance of the mpd code is very well pre- 
dicted by the model given in §2.3.4. 

The binned total astrometric errors are shown as black box- 
and- whiskers plots on the bottom panel of Fig.Q The total astro- 
metric error of an observation is the distance between the mea- 
sured (estimated) and true (actual) position of a star: Ar = 
■J (X — Xtmc) 2 + (y — 3^truo) 2 - The four grey limits seen on the 
bottom panel of Fig. Q are theoretical predictions (derived from 
§2.4.4) for the median (50% cumulative fraction: grey solid curve), 
top hinge (75%: bottom of the grey band), top fence (~98.97%: top 
of band), and 5-a outlier (99.99997%: grey dashed curve) values. 
The values of these theoretical limits are approximately equal to 
1.178 a y , 1.666 cry, 3.027c v, and 5.890 cry, where ay is the 
rms measurement error for the stellar X position. If the astrometric 
performance model is correct and mpd has been coded correctly, 
then (1) the observed median values should intersect the theoretical 
median value, (2) most of the top whiskers should be found inside 
the band, and (3) most of the outliers should be found above the 
top of the band and all of the outliers should found below the 5-cr 
outlier limit. 

Comparing the total astrometric errors of the 20,000 simulated 
CCD observations with the grey theoretical limits, one sees that the 
astrometric performance of the mpd code is very well predicted by 
the model given in §2.4.4. 

Figure|2|shows the relative stellar intensity errors and the rela- 
tive X position errors of the 20,000 stars analyzed in Fig.Q The rel- 
ative stellar intensity error is the difference between the measured 
(estimated) and true (actual) stellar intensity values divided by the 
estimated stellar intensity error: AS = (£ — £ true)/c£ . The rela- 
tive X position error is the difference between the measured (esti- 
mated) and true (actual) stellar X position values divided by the es- 
timated X error: AX = (X — A'truc) l a X- ^ m Pd has been coded 
correctly, the relative error distributions for the stellar parameters 
£ , X, and y should be normally distributed. The five grey limits 
seen on each panel are theoretical predictions (based on the normal 
distribution) for, from bottom to top, the bottom fence (~0.35% cu- 
mulative fraction: bottom of the bottom grey band), bottom hinge 
(25%: top of bottom band), median (50%: grey solid line at zero), 
top hinge (75%: bottom of top band), top fence (~99.65%: top 
of top band) values. If the relative errors for £ and X are indeed 
normally distributed, then (1) the observed median values should 
be near zero, (2) most of the whiskers should be found inside the 
bands, and (3) most of the outliers should be beyond the fence val- 
ues. 

Comparing the relative errors for £ and X of the 20,000 sim- 
ulated CCD observations with the grey theoretical limits, one sees 
that these errors are, as expected, normally distributed. 

• The mpd code works well with oversampled analytical Gaus- 
sian PSFs and its performance can be very well predicted with the 
photometric and astrometric models derived in §2. 



5.7.2 Discrete PSFs 

The 20,000 simulated CCD observations analyzed in Figs.Qand|2| 
were reanalyzed with mpd using an oversampled discrete Gaussian 
PSF with a FWHM of 3.0 px. Figure|3|shows the resultant absolute 
photometric errors and total astrometric errors. Figure|4]shows the 
resultant relative errors for £ and X. Notice how similar Figs.Qand 
|3|and Figs.|2|and|4]are to each other. 
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Figure 3. The absolute photometric errors (top) and total astrometric errors 
(bottom) of 20, 000 simulated CCD stellar observations used in Figs.Hland 
[2|were analyzed with mpd using an oversampled discrete Gaussian PSF 
with a FWHM of 3.0 px (/3 « 21.44 px 2 ; V = 1). 



Despite the very different way the shape information of the 
PSF was encoded (i.e., discrete versus analytical representations), 
mpd produced nearly identical photometric and astrometric results. 

How similar are the measured stellar positions? Figure [5] 
shows the relative X and 3^ position differences between the previ- 
ous analytical and numerical analyses with the mpd code. The top 
panel shows the difference between the numerical X result and the 
analytical X result divided by the estimated error of the analytical 
result. Similarly, the bottom panel shows the difference between 
the numerical 3^ result and the analytical 3^ result divided by the 
estimated error of the analytical result. The relative differences be- 
tween the numerical and analytical methods are not normally dis- 
tributed — observe how much smaller the values on the ordinate of 
Fig.|5|are compared to those of Figs.[2]and[4] Figs.[2|and|4]are nor- 
mally distributed and the source of the scatter is photon noise. Fig- 
ure [5] indicates that the relative differences between the numerical 
and analytical methods for astrometry are less than one-fifteenth 
the difference due to photon noise. In other words, the computa- 
tional noise due to the chosen analysis method (numerical versus 
analytical) is insignificant when compared to the unavoidable pho- 
ton noise due to the random arrival of photons in any astronomical 
CCD observation. 

• The mpd code works as well with oversampled discrete Gaus- 
sian PSFs as it does with oversampled analytical Gaussian PSFs. 
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Figure 4. Relative stellar intensity errors (top) and relative X position errors 
(bottom) of the data set used in Fig.l3l 
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Figure 5. Relative X and y position differences (top and bottom, respec- 
tively) between the numerical (subscript N) and analytical (subscript A) re- 
sults of the same 20,000 stars used in Figs. [T1l4l 



5.1.3 Inefficient Detectors 

While the volume, V, of the PRF was carefully tracked through- 
out the derivation of the photometric and astrometric performance 
models in §2, all previous simulations have assumed a perfect de- 
tector (V = 1). Let us now check to see if the effects of a PRF 
volume integral that is less than one has been correctly accounted 
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Figure 6. The absolute photometric errors (top) and total astrometric er- 
rors (bottom) of 20, 000 simulated CCD stellar observations analyzed with 
mpd using a discrete Gaussian PSF with a FWH M of 3.0 px with an inef- 
ficient detector with V = 1/9 (/3 w 1736.79 px 2 ). See the text for more 
details. 



for in the performance models of §2 by analyzing simulated obser- 
vations imaged on a very inefficient detector (V <C 1). 

Twenty thousand oversampled CCD stellar observations were 
simulated assuming a very inefficient detector with V = 1/9. Stars 
were simulated using a discrete Gaussian PSF with a FWHM = 3 
px located near the center of 60x60 pixels, the input stellar inten- 
sities ranged from —8 to —15 mag (1585 to 10 7); and the ob- 
served background sky level was assumed to be a constant value of 
B = 11.1111 e" (Btruc = 100 7, (V) = 1/9) all other simulation 
parameters were the same as before. 

All the simulated observations were analyzed with mpd in the 
same way as described for the numerical experiment shown in Fig. 
[5] — except that the volume of the PRF was set to V = 1/9 in order 
to simulate the use of an inefficient detector which converts only 
~ 1 1 . 1 % of photons to electrons. 

Figure[6|shows the absolute photometric errors and total astro- 
metric errors of this numerical experiment. The median effective- 
background area of PRF of these observations was (3 m 1736.79 
px 2 which is, as expected, 81 (= V - 2 ) times larger than the median 
value reported in Fig.0 

Comparing the simulation results with the grey theoretical 
limits, one sees that the photometric and astrometric performance 
of the mpd code is very well predicted by the theoretical perfor- 
mance models given in §2. 



The black dash-dot curves in each panel of Fig. |5| shows the 
expected median response with a perfect detector; these curves are 
the same as the solid grey median curves found in Fig. [3] The ob- 
served stellar intensities and observed background sky level are 
nine times fainter than was seen in the numerical experiment shown 
in Fig.|5|and the median photometric and astrometric errors in Fig. 
|6|are, as expected, ~3 (= V _1,/2 ) times larger when the inefficient 
detector is used. 



• The mpd code and the theoretical performance models work 
well with PRFs that have volumes of less than one. 



5.2 Undersampled Discrete PSFs 

Twenty thousand undersampled CCD stellar observations were 
simulated using an analytical Gaussian with a FWHM = 1.5 px, 
the other simulation parameters were the same as given in §5.1.1. 
The median effective-background area of PRF of these observa- 
tions was f3 « 6.12 px 2 (V = 1). All the simulated observa- 
tions were analyzed with mpd using a discrete Gaussian PSF with 
FWHM = 1.5 px. 

Figure^shows the absolute photometric errors and total astro- 
metric errors of this numerical experiment. While the photometric 
and astrometric results for stars with £true iS 30, 0000 e~ are fine, 
the results for stars brighter than that limit are seen to quickly de- 
grade in accuracy with the brightest stars having median errors that 
are ~40 times worse than expected. 

What starts going wrong at £tiue ~ 30,000 e~? Figure [8] 
shows a one pixel wide slice through a pixel-centered discrete 
Gaussian PSF with FWHM = 1.5 px that was shifted half of a 
pixel in X to the right using damped sine function given in equa- 
tion 1491 . The dashed black curve looks fine, but when expanded by 
a factor of 100, one sees that negative side lobes have been created 
due to the fact that the Nyquist-Shannon Sampling Theorem has 
been violated. Doing a sine interpolation (damped or otherwise) 
on undersampled data is never a good idea - the "ringing" seen in 
Fig.[3]is a classic signature of an edge that is too sharp to be ade- 
quately expressed with the limited spatial information contained in 
an undersampled observation. The biggest negative side lobe of the 
shifted PSF has a value of about —0.0006. Although that may seem 
to be a small value compared to the total volume integral of one, 
it is actually quite disastrous because negative PSF values have no 
physical meaning. 



It is now clear what has gone wrong for stars with £ true > 
30, 000 e~ . At stellar intensity values greater than 17,000 electrons, 
the intensity-scaled undersampled PSF models can have negative 
side lobes that are larger than the rms observed background sky 
noise level ( I -0.00061 * 17,000e" = 10.2 e" > 10 e~ » VB). 
At stellar intensity values greater than 167,000 electrons, the obser- 
vational models have physically nonsensical negative sky values. 

• Aliasing (ringing) effects will generally only be seen with 
bright stars since a large number of photons are required in order to 
adequately sample the higher spatial frequencies of a PSF. 

• Fitting undersampled observations of bright stars with under- 
sampled PSFs results in poor photometry and astrometry. 
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Figure 7. The absolute photometric errors (top) and total astrometric er- 
rors (bottom) of 20, 000 simulated CCD stellar observations analyzed with 
mpd using an undersampled discrete Gaussian PSF with a FWHM of 1.5 
px (P fa 6.12px 2 ;V = l). 
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Figure 9. The absolute photometric errors (top) and total astrometric er- 
rors (bottom) of 20, 000 simulated CCD stellar observations analyzed with 
mpd using a 2x2 supersampled discrete Gaussian PSF with a FWHM of 
1.5 px (/3 fa 6.17px 2 ;V = l). 




Figure 8. A one pixel wide slice through a pixel-centered discrete Gaussian 
PSF with FWHM = 1.5 px that was shifted half of a pixel in X to the right 
using equation The thick grey curve is the same PSF multiplied by a 
factor of 100. 



5.3 Supersampled Discrete PSFs 



A supersampled PSF is a PSF with pixels that have greater 
spatial resolution (higher spatial frequencies) than the actual pixels 
in the observational data. For example, a 2 x 2 supersampled PSF 
uses 4 pixels to describe every physical pixel of the CCD observa- 
tion; each supersampled pixel has twice the spatial resolution of the 
actual pixels in the observation. 

Twenty thousand undersampled CCD stellar observations 
were simulated using an analytical Gaussian with a FWHM = 1.5 
px; the other simulation parameters were the same as before. All 
the simulated observations were analyzed with mpd using a 2x2 
supersampled discrete Gaussian PSF with FWHM = 1.5 px (f3 m 
6.17 px 2 ; V= 1). 

Figure|9|shows the absolute photometric errors and total astro- 
metric errors of this numerical experiment. By providing mpd with 
extra information, in the form of a supersampled PSF, the Nyquist- 
Shannon Sampling Theorem was no longer violated and excellent 
photometry and astrometry was done with this undersampled data 
set. 

Figure [Tol shows the relative errors for £ and X. The rela- 
tive stellar intensity errors are normally distributed. However, the 



Stellar Photometry and Astrometry with Discrete PSFs 1 3 



4 

2 



I 
-2 
-4 



I 



x 



X 

I 

X 



4 
2 

-2 
-4 



'mm 



ni| — i i i ■ ■ i ■ ■ | . 

CM 



nr-i 



-H-H 



-H-H 



I I llllll I I M r 



TI 



t 



uL 



■J. 



■■I ' ■ 'it unit - 



10 s 10 3 10 4 10 6 10 6 
Electrons 

Figure 10. Relative stellar intensity errors (top) and relative X position 
eiTors (bottom) of the data set used in Fig.l9l 



relative X position errors are almost, but not quite, normally dis- 
tributed. The mpd code is accurately measuring the stellar posi- 
tions (i.e., the median difference, X — A'truc, values are zero) but 
the rms position error estimates (0%) are slightly underestimated 
(the top and bottom whiskers for Stmc > 10, 000 e~ are seen to 
extend beyond the grey bands). The same effect is seen with 3^- 
Using a higher-resolution supersampled PSF (3x3, 4x4, . . .) does 
not eliminate the small underestimation by mpd of position errors. 
The position errors estimated by mpd are close to the photonic 
limit, but the actual errors — for undersampled observa tions — 
are cl ose to the astrometric CRLB with square CCD pixels iWinicld 
1986). 



• Accurate and precise CCD stellar photometry and astrometry 
may be obtained with undersampled CCD observations if super- 
sampled PSFs are used during the PSF-fitting process. 



5.4 Critically-Sampled Discrete PSFs 

Let us now investigate what happens when critically-sampled data 
is fit with a critically-sampled PSF. 

Twenty thousand critically-sampled CCD stellar observations 
were simulated using an analytical Gaussian with a FWHM = 
2.35482 px; the other simulation parameters were the same as be- 
fore. All the simulated observations were analyzed with mpd us- 
ing a critically-sampled discrete Gaussian PSF with FWHM = 
2.35482 px 03 ^ 13.62 px 2 ; V = 1). 

Figure fTTI shows the absolute photometric errors and total as- 
trometric errors of this numerical experiment. Figure fT2l shows the 
relative errors for £ and X. Looking carefully at Figs. 1111 and 
1121 one sees that the photometric and astrometric performance is 
well matched to the theoretical expectations except for the bright- 
est three bins (£truc > 316, 000 e"). 
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Figure 11. The absolute photometric errors (top) and total astrometric er- 
rors (bottom) of 20, 000 simulated CCD stellar observations analyzed with 
mpd using critically-sampled discrete Gaussian PSF with a FWHM of 
2.35482 px (P Ri 13.62 px 2 ; V = 1). 
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Figure 12. Relative stellar intensity errors (top) and relative X position 
errors (bottom) of the data set used in Fig. II II 



Twenty thousand critically-sampled CCD stellar observations 
were simulated using an analytical Gaussian with a FWHM = 
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Figure 13. The absolute photometric errors (top) and total astrometric er- 
rors (bottom) of 20, 000 simulated CCD stellar observations analyzed with 
mpd using a 2x2 supersampled discrete Gaussian PSF with a FWHM of 
2.35482 px (P Ri 13.62 px 2 ;V = l). 



2.35482 px; the other simulation parameters were the same as be- 
fore. All the simulated observations were analyzed with mpd us- 
ing a 2x2 supersampled discrete Gaussian PSF with FWHM = 
2.35482 px (/3 w 13.62 px 2 ; V = 1). 

Fi sure [131 shows the absolute photometric errors and total as- 
trometric errors of this numerical experiment. Fieure fT4l shows the 
relative errors for£ and A 1 . The photometric and astrometric perfor- 
mance is well matched to the theoretical expectations for all stellar 
intensities. 



Comparing Fig. with Fig. [5] and Fig. II II with Fig. 1131 one 
sees that one can obtain excellent stellar photometry and astrometry 
with the MATPHOT algorithm for all stellar intensities - even if 
the observational data is undersampled - as long as the discrete 
PSFs used to do the model fitting are sampled finely enough to 
have sufficient spatial frequency coverage such that the Nyquist- 
Shannon Sampling Theorem is not violated. 

Comparing Fig.|5|with Fig. ll II one sees that the breakpoint for 
the MATPHOT algorithm between undersampled and oversampled 
data is 13.62 < f3 ^ 21.44 px 2 or, in terms of a Gaussian Full- 
Width-at Half Maximum, 2.35482 < FWHM sC 3 px. 
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Figure 14. Relative stellar intensity errors (top) and relative X position 
errors (bottom) of the data set used in Fig. ll 31 

one should use a supersampled discrete PSF which is oversam- 
pled in terms of supersampled pixels (spx). In other words, if the 
equivalent-background area is less than 21 square pixels ([3 < 21 
px 2 ; Gaussians: FWHM < 3.0 px), then one should use a super- 
sampled discrete PSF which has an equivalent-background area of 
at least 21 square supersampled pixels (f3 21 spx 2 ; Gaussians: 
FWHM > 3.0 spx). 

5.5 Ugly Discrete PSFs 

Let us now investigate the photometric and astrometric perfor- 
mance of the MATPHOT algorithm with an ugly (realistic) space- 
based PSF. 

Figure 1151 shows a simulated Next Generation Space Tele- 
scope (NGST) 1^-band CCD stellar observation. This simulated ob- 
servation used a 2x2 supersampled PSF which was based on a 8- 
meter TRW-concept 1.5 fi diffraction-limited primary mirror with 
1/13 wave rms errors at 1.5 fi; the original version of the PSF was 
kindly provided by John Krist (STScI). The six-sided "snowflake" 
pattern seen in Fig. l 1 5l is mainly due to fact that the primary mirror 
is composed of segmented hexagonal-shaped mirrors. Observers 
will note that this PSF is very similar to optical PSFs seen with the 
10-m telescopes at the W. M. Keck Observatory. The 6.5-m James 
Webb Space Telescope (JWST) is likely to have similar looking 
near-infrared PSFs once it achieves first light in ~201 1 . 

The NGST PSF is so complicated that it is unlikely that it 
could be represented adequately with a continuous analytical math- 
ematical function. Space-based observations frequently have high 
spatial frequencies which make them ideal candidates for photo- 
metric and astrometric analysis using discrete PSFs. 



Twenty thousand CCD stellar observations were simulated us- 
ing the simulated l/-band NGST 2x2 supersampled PSF described 
above; the other simulation parameters were the same as before. 
All the simulated observations were analyzed with mpd with the 
PSF used to create the simulated observations. 



If a discrete PSF is close to being critically sampled, then 



Photometry and Astrometry with Discrete PSFs 1 5 




: 




Figure 15. A simulated V-band Next Generation Space Telescope image 
based on a 2x2 supersampled PSF model for a 8-meter TRW-concept 1.5- 
micron diffraction-limited primary mirror with 1/13 rms wave errors. Con- 
tour levels of 90%, 50%, 10%, 1%, and 0.1% of the peak intensity are shown 
with black curves. The pixel scale is 0.0128 arcsec px _1 . This image uses 
a linear stretch with a pixel intensity mapping of black for < 70 e~ and 
white for > 150 e _ . 



Figure [Ttjl shows the absolute photometric errors and total as- 
trometric errors of this numerical experiment. Figure fTTl shows the 
relative errors for £ and X. The photometric and astrometric perfor- 
mance is well matched to the theoretical expectations for all stellar 
intensities. 

Although only Gaussian PSFs were used in previous numer- 
ical experiments, the excellent fit seen in the top panel of Fig. 1161 
between the theoretical photometric performance model (§2.3.4) 
and actual mpd measurements using such an ugly discrete PSF is 
not surprising once one remembers that the theoretical photometric 
performance model was derived from an abstract Point Response 
Function. 

The development of the theoretical astrometric performance 
model, however, required differentiation of the Point Response 
Function which I assumed to be an oversampled analytical Gaus- 
sian function. The analytical Gaussian bright star astrometric limit 
was transformed to the general form by assuming that the Gaussian- 
specific S 2 term could be replaced with the more general C 2 term, 
which, by definition, can be computed for any PRF. The same as- 
sumption was then used to derive the general faint star astrometric 
limit. The excellent fit seen in the bottom panel of Fig. ll6l indicates 
that this bold assumption is not only useful but practical. Many nu- 
merical experiments with very ugly discrete PSFs have shown that 
the theoretical astrometric performance model of §2.4.4 works well 
with ugly discrete Point Response Functions. 

If the MATPHOT algorithm is optimally extracting photomet- 
ric and astrometric information from a stellar observation, and mpd 
has been correctly coded, and the CCD observation has been prop- 
erly calibrated, and the PRF used in the observational model is cor- 
rect, and accurate estimates of the measurements errors for each 
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Figure 16. The absolute photometric errors (top) and total astrometric 
errors (bottom) of 20, 000 simulated CCD stellar observations analyzed 
with mpd using the 2x2 supersampled NGST PSF described in Fig.[T5l 
(J3 a 31.25 px 2 ;V=l). 
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Figure 17. Relative stellar intensity errors (top) and relative X position 
errors (bottom) of the data set used in Fig. 1161 
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Figure 18. A comparison between the cumulative \ 2 distribution for 3596 
degrees of freedom (thick curve) and the measured \ 2 value (thin curve), 
of the data set used in Fig. 1 161 reported by the mpd implementation of the 
MATPHOT algorithm. 
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Figure 19. The measured electron loss of the 10,000 simulated CCD obser- 
vations of —13 mag stars analyzed with mpdx using a 4x4 supersampled 
version of the NGST PSF. The electron loss is plotted as a function of the 
absolute value of the distance from the center of a star and the center of the 
active region of the central pixel of the stellar image. 



pixel have been made, then one expects that the \ 2 goodness-of-fit 
value reported by mpd to be distributed as a \ 2 distribution with 
the number of degrees of freedom equal to the difference between 
the number of pixels in the observation and the number of free pa- 
rameters. Figure fT8l shows that this prediction about the precision 
and accuracy of the MATPHOT algorithm has been verified: the cu- 
mulative distribution of the \ 2 reported by mpd (thin black curve) 
is seen to lie on top of the cumulative \ 2 distribution of for 3596 
[= 60 2 pixels — 4 free parameters {£, X, y, and B) ] degrees of 
freedom (thick grey curve). 



• The mpd code works well with ugly discrete PSFs and its per- 
formance can be well predicted using the general theoretical pho- 
tometric and astrometric performance models given in §2. 

• The x 2 goodness-of-fit value reported by mpd is a statisti- 
cally reliable measure of the quality of a photometric and astromet- 
ric reduction of a stellar observation obtained with the MATPHOT 
algorithm using ugly discrete PSFs. 

5.6 Ugly Detectors 

Let us now investigate the photometric and astrometric perfor- 
mance of the MATPHOT algorithm with an ugly PSF and an ugly 
detector. 

Suppose one has a detector where every pixel has been divided 
up into 16 square regions called "subpixels". Let us call the first row 
and first column of subpixels "gate structures" which are optically 
inactive with 0% QE. The remaining 9 subpixels are the optically 
active part of the pixel with a 100% QE. By definition, such a pixel 
would have a very large intrapixel QE variation with only 56.25% 
of the total pixel area being capable of converting photons to elec- 
trons. 

A few extra lines of code were added to the mpd program to 
simulate the image formation process with such an ugly detector. 
The new version of mpd is called mpdx and was designed specif- 
ically for use with 4x4 supersampled PSFs. 



Ten thousand CCD stellar observations of —13 mag stars 
(~2.512 13 7) were simulated and analyzed with mpdx using a 
4x4 supersampled version of the simulated IZ-band NGST PSF 
described above. The observed background level was assumed to 
be a constant value of B = 56.25 e~ (,6truc = 100 7, (V) = 
0.5625) and all other simulation parameters were the same as be- 
fore. The measured PRF volume of these simulated observations 
was V = 0.5616 ± 0.0185 which is consistent with the expected 
value of 0.5625 from the physical structure of a single pixel. The 
median and semiquartile range of the effective-background area 
(/3) of these observations was, respectively, 28.10 and 4.82 px 2 . 
The median critical-sampling scale length of these observations 
was C ~ 0.8398 px — indicating that these observations were 
undersampled, as expected. 

The optically inactive gate structures of the pixel cause the ob- 
served number of electrons in each stellar image to be significantly 
less than the number of photons which fell on the detector. The to- 
tal amount of loss was dependent on where the center of the star 
fell within the central pixel of the stellar image. Figure fl9l shows 
that stars centered in the middle of the active area of a pixel suf- 
fered a ~40% loss (Am ~ 0.56 mag) while those centered on gate 
structures (grey points) lost up to 47% (Am ~ 0.69 mag). 

Although this numerical experiment may seem to be very arti- 
ficial, large intrapixel sensitivity variations can be found i n cam- 
eras c urrently installed on the Hubble Space Telescope. iLaued 
1 1999) reported peak-to-peak variations of 0.39 mag at the J band 
(Fl 10W) and 0.22 mag at the H band (F160W) of the NIC3 cam- 
era of the HSTNICMOS instrument. The peak-to-peak variation of 
^0.2 mag at F160W wi th NIC 3 was independently confirmed by 
iHook & FruchteiH2000l) . 



The mean observed stellar magnitude for these —13 mag stars 
was — 12.3728±0.0359 mag. The photometric performance model 
predicts an rms measurement error of 0.0036 mag for these bright 
stars. With an average loss of 44% and an rms measurement error 
that is ten times larger than expected from photon statistics, the 
observed stellar magnitudes were neither precise or accurate. 
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Figure 20. The observed (left) and the measured (right) stellar magnitude 
distributions of the 10,000 —13 mag stars described in Fig. 1191 



Figure l20l shows that mpdx was able to do an excellent job in 
recovering the true stellar magnitude of the 10,000 —13 mag stars 
— despite being presented with a worst-case scenario of under- 
sampled observations with an ugly PSF imaged on an ugly detector 
with a very large intrapixel QE variation. 



The mean measured stellar magnitude reported by mpdx was 
— 12.9998 ± 0.0039 mag and the mean rms error estimated by 
mpdx was 0.00384± 0.00006 mag. The photometric performance 
of mpdx is fully consistent with theoretical expectations — which 
were derived for an ideal detector with no intrapixel QE variation. 

Twenty thousand CCD stellar observations were simulated 
and analyzed with mpdx using the same 4x4 supersampled ver- 
sion of the simulated V-band NGST PSF. The input stellar inten- 
sities ranged from —6 to —15 mag (251 to 10 6 7). The observed 
background level was assumed to be a constant value of B — 56.25 
e~ (Btrue = 100 7, (V) = 0.5625) and all other simulation pa- 
rameters were the same as before. The median and semiquartile 
range of the effective-background area (/3) of these observations 
was, respectively, 28.04 and 4.77 px 2 . 

Figure ITTl shows the absolute photometric errors and total as- 
trometric errors of this numerical experiment. Comparing the simu- 
lation results with the grey theoretical limits, one sees that the pho- 
tometric and astrometric performance of the mpdx code is well 
predicted by the theoretical performance models given in §2. 



• Excellent stellar photometry and astrometry is possible with 
ugly PSFs imaged onto ugly detectors as long as the image for- 
mation process within the detector is accurately modeled by the 
photometric reduction code. 



6 DISCUSSION 

After developing theoretical photometric and astrometric perfor- 
mance model for Point Spread Function (PSF)-fitting stellar pho- 
tometry, I described the unique features of the MATPHOT algo- 
rithm for accurate and precise stellar photometry and astrometry us- 
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Figure 21. The absolute photometric errors (top) and total astrometric er- 
rors (bottom) of 20, 000 simulated CCD stellar observations analyzed with 
mpdx using a 4x4 supersampled version of the NGST PSF (f3 ^ 28.04 
px 2 ;V = 0.5625). 



ing discrete Point Spread Functions. I conducted numerical exper- 
iments with the mpd implementation of the MATPHOT algorithm 
and demonstrated that the computational noise due to the chosen 
analysis method (numerical versus analytical) is insignificant when 
compared to the unavoidable photon noise due to the random arrival 
photons in any astronomical CCD observation. The MATPHOT al- 
gorithm was specifically designed for use with space-based stellar 
observations where PSFs of space-based cameras frequently have 
significant amounts of power at higher spatial frequencies. Using 
simulated NGST CCD observations, I demonstrated millipixel rel- 
ative astrometry and millimag photometry is possible with very 
complicated space-based discrete PSFs. 

The careful reader will observe that I have not discussed how 
a discrete PSF is derived. The MATPHOT algorithm will optimally 
determine the brightness and position of a star in a CCD obser- 
vation when provided with the correct PSF and DRF — functions 
which need to be determined beforehand through calibration pro- 
cedures. Photometric and astrometric accuracy and precision de- 
grades if either the PSF or DRF is poorly known. PSF reconstruc- 
tion (calibration) is a complicated topic in its own right and has 
been the subject of many articles, instrumentation reports, and en- 
tire workshops. The challenges of PSF reconstruction are many. An 
astronomer may be faced with trying to derive a PSF from an ob- 
servation 
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• with a variable PSF within the field of view, 

• that has too few bright stars, 

• that might be undersampled, 

• that might be poorly dithered, 

• that might be poorly fiat-fielded, 

• that exhibits significant charge transfer efficiency variations, 

• that has variable charge diffusion within the CCD substrate, 

• with significant photon loss due to charge leakage, 

• that might not actually be linear below the 1% level. 

While many of these problems can be overcome by the proper de- 
sign of instruments or experiments, their solution is beyond the 
scope of this article which has sought to determine the practical 
limits of PSF-fitting stellar photometry. 

The analysis presented in this article has assumed that PSFs 
are perfectly known - a situation that is rarely, if ever, physically 
possible. The cores of observationally based PSFs are generally 
much better determined than the broad wings due to simple pho- 
ton statistics. The effect of large instrumental calibration errors can 
also be significant. For example, flat-field limitations can dramat- 
ically impact the achievable levels of photometric and astrometric 
precision. An investigation based on theory of the effect of PSF er- 
rors and flat-field calibration error on the limits of PSF-fitting stel- 
lar photometry would be very difficult. An investigation based on 
numerical experiments, however, might be a much more tractable 
proposition. In any case, a through investigation of the effects of 
calibration errors on the limits of PSF-fitting stellar photometry is 
best left to another article. 

I would like to thank the anonymous referee whose careful 
reading and insightful comments have improved this article. Spe- 
cial thanks are due to Ian Roederer who checked all of the math- 
ematical proofs in an early draft of the manuscript. Jessica Moy 
has my heartfelt thanks for cheerfully helping me acquire copies 
of many articles in technical journals that are not available in the 
NOAO library collection. I would also like to thank the following 
people for providing stimulating conversations, encouragement, 
and support during the extensive research and development effort 
behind the creation of the MATPHOT algorithm and and its various 
software implementations: Joseph Bredekamp, Susan Hoban, Dot 
Appleman, Taft Armandroff, Nick Buchholz, Marc Buie, Harvey 
Butcher, Julian Christou, Chuck Claver, Lindsey Davis, Michele De 
Le Pena, Mike Fitzpatrick, Ken Freeman, Dan Golombek, Richard 
Green, Steve Howell, George Jacoby, Buell Jannuzi, Stuart Jef- 
feries, Ivan King, John Krist, Todd Lauer, John MacKenty, John 
Mather, Mike Merrill, Dave Monet, Mark Morris, Jeremy Mould, 
Jan Noordam, John Norris, Sudhakar Prasad, the late Alex Rogers, 
Steve Ridgway, Pat Seitzer, James Schombert, Donald West, and 
Sidney Wolff. This research has been supported in part by the fol- 
lowing organizations, institutions, and research grants (in reverse 
chronological order): National Aeronautics and Space Adminis- 
tration (NASA), Interagency Order No. S-13811-G, which was 
awarded by the Applied Information Systems Research Program 
(NRA 01-OSS-01) of NASA's Science Mission Directorate, Inter- 
agency Order No. S-67046-F which was awarded the Long-Term 
Space Astrophysics Program (NRA 95-OSS-16) of NASA's Office 
of Space Science, Center for Adaptive Optics (CfAO), Institute for 
Pure and Applied Mathematics (IPAM), Kitt Peak National Ob- 
servatory, National Optical Astronomy Observatory, Association 
of Universities for Research in Astronomy Inc. (AURA), National 
Science Foundation, Mount Stromlo and Siding Spring Observa- 
tories, Australian National University, the Netherlands Foundation 



for Astronomical Research (ASTRON), and the Netherlands Orga- 
nization for the Advancement of Pure Research (ZWO), and the 
Kapteyn Astronomical Institute of Groningen. 



REFERENCES 

Abramowitz, M., & Stegun, I. 1964, Handbook of Mathematical 
Functions with Formulas, Graphs, and Mathematical Tables, Ap- 
plied Mathematics Series, 55, eds. M. Abromowitz and I. Stegun 
(Washington, D.C.: NBS) 

Arfken, G. B. 1970, Mathematical Methods for Physicists, 2nd 
ed., (New York: Academic Press) 

Biretta, J., et al. 2001, WFPC2 Instrument Handbook, Version 6.0 
(Baltimore: STScI) 

Hook, R. N., & Fruchter, A. S. 2000, in ASP Conf. Ser., Vol. 216, 
Astronomical Data Analysis Software and Systems IX, eds. N. 
Manset, C. Veillet, and D. Crabtree (San Francisco: ASP), 521 

Irwin, M. J. 1985, MNRAS, 214, 575 

Jakobsen, P., Greenfield, P., & Jedrzejewski, R. 1992, A&A, 253, 
329 

King, I.R. 1971, PASP, 83, 199 
King, I.R. 1983, PASP, 95, 163 
Lauer, T.R. 1999, PASP, 111, 1434 
Levenberg, K. 1944, Q. App. Math., 2, 164 
Marquardt, D. 1963, J. SIAM, 11, 431 
Mighell, K. J. 1989, MNRAS, 238, 807 

Mighell, K. J. 1999, in ASP Conf. Ser., Vol. 172, Astronomical 
Data Analysis Software and Systems VIII, eds.D. M. Mehringer, 
R. L. Plante, and D. A. Roberts (San Francisco: ASP), 317 

Mighell, K. J. 2002, in Proceedings of SPIE, Vol. 4847, Astro- 
nomical Data Analysis II, eds. J.-L. Starck and F. D. Murtagh 
(Bellingham, WA: SPIE), 207 

Pogson, N. 1856, MNRAS, 17, 12 

Perryman, M. A. C, Jakobsen, P., Colina, L., Lelievre, G., Mac- 
chetto, F, Nieto, J. L., & di Serego Alighieri, S. 1989, A&A, 
215, 195 

Press, W. H., Flannery, B. P., Teukolsky, S. A., & Vetterrling 

W. T. 1986, Numerical Recipes. The Art of Scientific Computing 

(Cambridge: Cambridge University Press) 
Tukey, J. W. 1977, Exploratory Data Analysis, (Reading, MA: 

Addison- Wesley) 
Wells, D. C, Greisen, E. W., & Harten, R. H. 1981, A&AS, 44, 

363 

Winick, K. A. 1986, JOS A A, 3, 1809 



APPENDIX A: BOX- AND- WHISKER PLOTS 

A box-and-whisker plot (a.k.a. box plot) is a graphical method of 
showing a data distribution. A box is drawn showing the inner quar- 
tile range of the data which, by definition, includes half of all the 
data values (see Fig. lAH . The median of the data is shown with a 
bar inside the box. The bottom end of the box is the lower quar- 
tile (25%) of the data: iTukevI Jl977l) . the creator of the box-and- 
whiskers plot, calls this value the lower hinge, LH, value. The top 
end of the box is the upper quartile (75%) of the data which is 
called the upper hinge, UH, value. The step value is 1.5 times the 
inner quartile range: A = 1.5 * (UH — LH). The top fence value 
is the sum of the upper hinge and step values: TF = UH + A. The 
bottom fence value is the difference between the lower hinge and 
step values: BF = LH — A. The top whisker is drawn from the 
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Figure Al. A box-and- whiskers plot of a data set of 1000 normal deviates. 
See the text for details. 



upper hinge value to the largest data value that is less than or equal 
to the top fence value: TW TF. Similarly, the bottom -whisker is 
drawn from the lower hinge value to the smallest data value that is 
greater than or equal to the bottom fence value: BW ^ BF. Data 
values that are greater than the top fence value or less than the bot- 
tom fence value are called outliers and are plotted at their appro- 
priate value beyond the whiskers. For a normal distribution, which 
is a Gaussian distribution with a mean of zero and a standard de- 
viation of one, the bottom fence, bottom hinge, median, top hinge 
and top fence values are, respectively, —2.6980 (0.35% cumula- 
tive fraction), -0.6745 (25%), zero (50%), 0.6745 (75%), 2.6980 
(99.65%). 



Figure lATl shows a data set of 1000 normal deviates. The his- 
togram of the data with 0.1-wide bins is shown with thin black 
lines. The cumulative fraction distribution of the data is shown as 
a thick gray curve. The box-and- whisker plot of the data is shown 
with thick black lines below the histogram; arrows show the re- 
lationship between various box values and the cumulative frac- 
tion distribution. The mean and standard deviation of this data set 
are, —0.0341 and 0.9739, respectively. The bottom fence, bottom 
whisker, bottom hinge, median, top hinge, top whisker, and top 
fence values of this data set are, respectively, —2.5511 (0.25% cu- 
mulative fraction), -2.4940 (0.30%), -0.6522 (25.10%), -0.0231 
(50.00%), 0.6137 (75.10%), 2.4580 (99.50%), 2.5126 (99.57%). 
The seven outlier values of this data set, —2.9500, —2.6320, 
2.5390, 2.7150, 2.7430, 2.8270, 2.8530, are plotted in Fig.lATlas 
diamonds beyond the whiskers. 



